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ABSTRACT 


Numerical solution of fractional Lane-Emden equation which is 
the nonlinear singular initial value problem in astrophysics has 
been proposed by fractional order of Rationalized Haar operational 
matrices based on collocation method. The advantage of our 
technique is that the computational speed is high due to using the 
RH operational matrices and the convergence rate is exponential. 
We investigate two examples to test the accuracy, speed and 
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INTRODUCTION 


Fractional differential equations are generalizations 
of classical differential equations, which are 
obtained by replacing integer-order derivatives by 
fractional ones. Their advantages comparing with 
integer-order partial differential equations are the 
capability of simulating natural physical process 
and dynamic system more accurately. Furthermore, 
the fractional order models of real systems are 
regularly more adequate than usually used integer 
order models [1]. 


The ordinary Lane-Emden equation describes the 
temperature variation of a spherical gas cloud under 
the mutual attraction of its molecules and published 
by Lane [2] and further explored in detail by 
Emden [3]. 


Lane-Emden equation is a second-order ordinary 
differential equation with an arbitrary index, known 
as the polytropic index, that deals with the issue of 
energy transport, through the transfer of material 
between different levels of the star [4]. 


The ordinary Lane-Emden does not describe the 
dynamics of systems in complex media accurately. 
Also, the memory effects are better described 
within the fractional derivatives, therefore the 
fractional Lane-Emden equations extract hidden 
aspects of the complex phenomena and it can 
describe the parameters better and more accurately 


[3]. 


In this paper, we study the Caputo fractional 
derivative of Lane-Emden type equation as [6]: 


a 


Déy(t)+z—DiyH+fiy=90, 0< 
t<b, Aaz=0, (1) 

with initial conditions 

y(0)=A, y'(0)=B, (2) 


where 0<a<2,0<PB <1 and p#, DP refer to 
the Caputo fractional derivative of order a, 6 with 
respect to t, and parameters A,B and A is constant, 
f(t, y) is a nonlinear function of t, y and g(t) is a 
continuous real valued function in C[0, b). 
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The single standard and general solution algorithm 
for problems regarding fractional calculus has not 
yet been constituted and finding affirmative, 
accurate, stable and credible solution methods along 
with fast time computation order is beneficial. 


For the solution of nonlinear and singular fractional 
Lane-Emden equations which are difficult to solve 
there are many numerical and analytical methods 
such as, differential transform method [6], 
collocation technique [7], finite difference 
algorithm[8], Homotopy analysis approximation [9], 
reproducing kernel method [10], hybrid wavelet 
method [11], Legendre wavelets spectral technique 
[12], coupling of wavelets and Laplace transform 
[13], spectral Legendre’s derivative algorithms [14], 
discontinuous finite element approximation [15] and 
Galerkin operational matrices [16]. 


In this paper, RH collocation method based on 
fractional order of operational matrices is developed 
to obtain solutions of singular fractional Lane— 
Emden type equations. The properties of the RH 
orthogonal functions are used to convert the 
problem into a system of algebraic equations which 
can be solved by suitable algorithms with 
exponential convergence rate. 


The organization of the rest of this paper is as 
follows: In Section 2, we describe a_ short 
introduction to the basic definitions of the fractional 
calculus and Rationalized Haar properties which 
drive some tools for developing our method. In 
Section 3, we summarize the application of the RH 
operational matrices collocation method for solving 
the model equation. In Section 4, the proposed 
method is applied to some types of fractional Lane- 
Emden equations, and comparisons are made with 
the existing analytic solutions that were reported in 
other published works in the literature. The 
conclusions are described in the final section. 


Basic definitions 

A. Definitions of fractional derivatives and 
integrals 

In this section, we start with recalling the essentials 

of the fractional calculus [1]. 


Definition 1. For m to be the smallest integer that 
exceeds a, Caputo’s time-fractional derivative 
operator of order a > 0 is defined as: 

pre-e D m y(t) 
Cha t t , 

bee =| 

a’t Y ( ) pr y(t), 
where 


Ify(@) = eG —t)*1y(t)dt, t,a>0,(4) 


m-1<a<m, 
a=mmeNn, 


(3) 


Some of the most important properties of operator 
If, Df is as follows: 


are) — Putt) ita 
T(t") aT eae 


PO) =a Te"); 


T(ut1) a 
Ad aire all 


Di (xy) + ps(t)) = KDE y(t) + pDEs(), 
where xk and p are constants. 


J 


(5) 


J 


B. Rationalized Haar functions 
The orthogonal set of RH functions is a group of 
square waves with magnitude of +1 in some 
intervals and zeros elsewhere [17]. The RH 
functions is defined on the interval C = [0, b) by: 
2 

A= 4 Lbst< ly (6) 

2 
0, otherwise 
where 

se = 9:4 
lg = tb, q=0,5,1.(7) 


The value of u is defined by two parameters m and 
Nn as: 


u=2™+n—-1, m=0,1,2,..., n=1,2,3,...,2™. 
(8) 

ho(t) is defined form = n = 0 and is given by: 
ho(t)=0, O<St<b.(9) 


We can expand any function y(t) € L?[0, b) in first 
k terms of RH functions as: 


y(t) = YE xy hy (t) =X"; (0), (10) 
where 
c= 2s =012 x 


The RH functions coefficient vector X and RH 
functions vector ¢;(t) are defined as: 
X=[Xo, Xy, wy Xp-1]", (11) 


b(t) = [ho(t), hi(t), «, Ayes (0)]", 12) 


The matrix ®;,., can be expressed as: 


1 3 2k-1 
Poe = | bu (Fe) Px (Se) be (GE) 9) 
C. Rationalized Haar operational matrices of 
fractional order integration 
The integration of the @,(t) defined in Eq. (12) is 
given by: 


fo bx (s)ds ~ Jp, (t), (14) 


where J = Jxy% is the k X k operational matrix for 
integration and is given in [18] as: 
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| “Oxo OO 


leq s, on 0 , (15) 
(2)x(@) 
where ®, = [1], h= [5]. 


Also, we obtain the Rationalized Haar operational 


matrix of the fractional order integration 
Ci bx) (t) [19] as: 
CEO O) = Ske Px (O), (16) 
where 
Tix = Perk *Pyexe, (17) 
where 
1 Wy, @2 Wr-1 
. 4 O 1 @y . Wr-2 
a —_-" 
7 ee Grea 0 0 1 ..) WR-3} 8) 
0 O * ; 
0 0 O O 1 
With @, = (r+ 1)8t? —2re** + @ — 18 


Theorem. 1 Assume that y(t) is continuous and 
bounded on (0,1), then there is [20] : 


ID>0, wee (01) |y™ | < D. (19) 


Then, the error norm at Tth level of resolution 
satisfies the following inequality: 


I De yt) — De yr(t) les 
TP-2 P(p— i (p- a) [1=22(@-p) a p) oi 


where D? y(t) is the following approximation of 
Dfy(), 
Di yr(t) = 


and 


y(t) lle= (fp y? (de)*. 22) 


Implementation of Rationalize Haar operational 
matrices for solving the model 

To solve Eq. (1) with boundary conditions in Eq. 
(2) we let: 


y"(t) = Lez Xu hu (t) = Xb, (0), (23) 
where 
u=2™+n-1, 
(24) 

Using Eqs. (14) and (23) we get: 

y(t) = X"Jo,() + B, (25) 

and 

y(t) = XJ fo by (de + Bt + A = XTJ?p,(t) + 
Bt + A. (26) 


u=0 Xu Ay (E), (21) 


m=0,1,2,..., n=1,2,3,...,2 


Also by using Eq.(5), for 1 <a <2,0< 6 <1, we 
a B,_ ace a B 

DFt = 0,D,t = FO=py! ,DfA = D, A=), 

and using Eqs.(26), (3) we have: 

Df y(t) = X7J?2-% db, (t) + BDEt + DEA = 

XTC“ ,(€), (27) 

De y(t) = OP Gel +BDet+DeA= 

XTJ2-P) gb, (t) + ——. (28) 


have 


ce B) 
Using Eqs. (1), (27) and (28) we get: 


Xe Dh, (t) ak = “5 (x72 (0) + 
aan —) + f(t, X7J*,(t) + Bt + A) = g(t). (29) 


The residual Res; (t) for Eq. (1) can be written as: 
Res;,(t) = 


t1-B 
i 2= ow) 
+f (t,X7J7o,(t) + Bt + A) — g(t). (30) 


The equations for obtaining the coefficients x, 
arises from equalizing Res,(t) to zero at k RH 
collocation points defined by: 


t= bp, 1=1,2,...,k. (31) 


2K 


A (x7 6.(0 + 72 + 


By substitution collocation points t;,i = 1,2,...,k 
in Res;,(t) and equalizing to zero we have: 


Res,(t;) =0, i=1,2,...,k. (32) 


Eq. (32) gives k nonlinear algebraic equations 
which can be solved for the unknown coefficients 
xX,u=0,1,..,k—-—1 by using the well-known 
Newton’s method. Consequently, y(t) given in Eq. 
(1) can be calculated. 


Numerical results and illustrative examples 

In this section, we solve the model numerically by 
the proposed method and show the efficiency of the 
method with the numerical results of two examples. 


Example 1. We consider the following fractional 
equation which has been solved in [6],[7]: 


m_ Déy(t) + Dey(t) + zz y(t) = s(t), (33) 


with initial conditions 

y(0) =0, y'(0) = 0, (34) 
where s(t) =t?~? (6c (“+ 
2 + Tearew 

The exact solution to this problem for a = 1.5, 8 = 
0.5 andd = 1ist? —t?. 


4 LG=B)+Ar(4= ©) 
r(4-B)r(4-a) 
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Table 1 shows the comparison of the y(t) between 
numerical solutions obtained by the method 
proposed in this letter for k = 8,16,32,64 and their 
absolute errors with respect to the exact solution. 
Fig. 1. shows the numerical and analytic graph of 
solution of Eq.(33 ) for k = 32. 
The graph of absolute errors for k = 32 at RH 
i-1 


. ; 2 : . 
collocation points t; = je 1,2,...,k is shown 


in Fig. 2. 
Example 2. Consider to the following equation [6], 
[7]: 
a 1 
Di y(t) + sap DEVO) + aa yD = rO), BS) 


with initial conditions 
y(0) =0, y'() =0, (36) 


_ realy (2 4 LG=B)+ar-a)) _ 

where r(t)=t (2 (- + 7G-OrG-a) ) 
6 (5 (4-B)+ar (4-a@) 

6 | (4—B)r(4—a) 
this problem for a=1.5,8 =1 and A=1 is 
t?7—t?. 
In table 2, we report the absolute errors with respect 
to the exact solution for different values k. 
In Fig. 3, the graph of y(t) for k = 32 and exact 
solution is plotted. 


)) The exact solution to 


The L*—Error and L®—Error are used to 
explore the dependence of errors on the parameters 
k. Fig. 4 displays the logarithmic scale of error 


versus k for k = 4,8,16,32,64,128 at RH 
collocation points. Since in this  semi-log 
representation the error variations are 
Conclusion 


approximately linear versus k, we observe that the 
values of error decay exponentially. 


Example 3.Consider to the following equation 
[21]: 


D2%y(t) + * Dey(t) +1=0,0 <a <1, (37) 


with initial conditions 

y(0) =1, y'(O) = 0. (38) 

The exact solution to this problem is 1— =. 
Table 3 displays the max absolute errors with 
respect to the exact solution for different values of 
k,a at RH collocation points. 
We can see clearly that better accuracy can be 
achieved by increasing the values of k and the 
values of error decrease when the values of a@ 
converge to 1. Fig. 5 shows the approximate 
solution obtained by present method for k = 64 for 
different values of a. 


Example 4.Consider to the following equation 
[21]: 


D?*y(t) + = Déy(t) + e”(t) = 0,0 < a < 1, (39) 


with initial conditions 
y(0)=1, y'(O) =0. (40) 


: ; . t2% 
The exact solution to this problem is = + 
t+e pou 


120a* 189006" 

Fig. 6 displays the behavior of the solution at 
a = 0.5,0.75,0.85,1 for k = 64. Fig. 7 displays the 
logarithmic scale of error versus k for k = 


8,16,32,64 at RH collocation points. 


RH operational matrices have been successfully employed to obtain accurate numerical solutions of 
fractional Lane—Emden type equations with linear and nonlinear terms, with fast convergence rate. The 
comparison between numerical and analytic methods has been made and it can be clearly seen from absolute 
errors in tables 1,2 and 3 that the numerical results of the suggested method converges to the exact solution 
by increasing the values of k. In addition, the properties of the RH functions have been used to convert the 
problem into a system of algebraic equations which can be conveniently solved by suitable algorithms. 


Table 1 Absolute errors with respect to the exact solution for k = 8, 16,32, 64 for Example 1 


k\t 0.2 0.4 0.6 0.8 1 
8 | 1.165e-3 | 1.863e-3 | 1.793e-3 | 1.539e-3 | 1.161e-3 
16 | 2.105e-4 | 2.256e-4 | 2.098e-4 | 1.759e-4 | 1.303e-4 
32 | 1.017e-6 | 1.377e-6 | 1.437e-6 | 1.292e-6 | 1.004e-6 
64 | 1.872e-8 | 1.997e-8 | 1.840e-8 | 1.526e-8 | 1.118e-8 
Table 2. Absolute errors with respect to the exact solution for k = 8,16, 32,64 for Example 2. 
k\t 0.2 0.4 0.6 0.8 1 
8 | 1.719e-3 | 1.924e-3 | 1.840e-3 | 1.572e-3 | 1.182e-3 
16 | 2.174e-4 | 2.317e-4 | 2.145e-4 | 1.794e-4 | 1.326e-4 
32 | 1.110e-6 | 1.443e-6 | 1.483e-6 | 1.324e-6 | 1.024e-6 
64 | 1.922e-8 | 2.038e-8 | 1.869e-8 | 1.547e-8 | 1.131e-8 
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Table 3. Max absolute errors with respect to the exact solution for k = 8, 16, 32, 64 for different 


values of a at RH collocation points for Example 3. 


k\a 0.25 0.5 0.75 0.85 1 
8 2.004e-3 | 1.735e-3 | 1.254e-3 | 1.193e-3 | 1.126e-3 
16 | 4.019e-4 | 3.758e-4 | 3.024e-4 | 1.994e-4 | 1.728e-4 
32 | 3.128e-6 | 2.854e-6 | 2.019e-6 | 1.917e-6 | 1.520e-6 
64 | 2.954e-8 | 2.158e-8 | 1.978e-8 | 1.287e-8 | 1.015e-8 
0 
-0.1 
yt) 
-0.2 
-0.3 
0 0.2 0.4 0.6 0.8 l 


‘ 


Fig. 1. Numerical (solid circle) and analytic (solid line) graph of solution of Example. 1, for k = 32. 
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1.2 =x 107° 


8.x 107! 
6.x 1077 
4.x 107" 
2.x 1077 
0 0.2 0.4 0.6 0.8 1 
t 
Fig. 2. Absolute errors at RH collocation points t; = —, i= 1,2,...,k, for Example. 1, for k = 32. 


@ IJTSRD | Unique Paper ID — NTSRD52389 | Volume -—6 | Issue—7 | November-December 2022 Page 700 


International Journal of Trend in Scientific Research and Development @ www.ijtsrd.com eISSN: 2456-6470 


0.3 


ya) 


0.1 


0 0.2 0.4 0.6 0.8 l 


t 
Fig. 3. Numerical (solid circle) and analytic (solid line) graph of solution of Example. 2, for K=32. 
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Fig. 4. Plot of the logarithmic scale of L? and L®- Error versus k, for Example 2. 
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Fig. 5. Approximation of the solution for different values of a for k = 64 for Example 3. 
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O8 


Fig. 6. Approximation of the solution for different values of a for k = 64 for Example 4. 
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Fig. 7. Plot of the logarithmic scale of L? and L”- Error versus k, for for Example 4. 
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